Phase field theory of crystal nucleation in hard sphere liquid 
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The phase field theory of crystal nucleation described in [L. Granasy, T. Borzsonyi, T. Pusztai, Phys. Rev. 
Lett. 88, 206105 (2002)] is applied for nucleation in hard-sphere liquids. The exact thermodynamics from 
molecular dynamics is used. The interface thickness for phase field is evaluated from the cross-interfacial 
variation of the height of the singlet density peaks. The model parameters are fixed in equilibrium so that the 
free energy and thickness of the (111), (110), and (100) interfaces from molecular dynamics are recovered. 
The density profiles predicted without adjustable parameters are in a good agreement with the filtered densities 
from the simulations. Assuming spherical symmetry, we evaluate the height of the nucleation barrier and the 
Tolman length without adjustable parameters. The barrier heights calculated with the properties of the (111) and 
(110) interfaces envelope the Monte Carlo results, while those obtained with the average interface properties fall 
very close to the exact values. In contrast, the classical sharp interface model considerably underestimates the 
height of the nucleation barrier. We find that the Tolman length is positive for small clusters and decreases with 
increasing size, a trend consistent with computer simulations. 



I. INTRODUCTION 

In this paper we present a detailed quantitative test of the 
phase field theory of crystal nucleation in the hard sphere (HS) 
fluid, a case in which all the data necessary for such analysis 
are available. 

Probably the least understood stage of the freezing of liq- 
uids is the initial nucleation phase in which the crystalline 
phase appears in the homogeneous liquid via heterophase fluc- 
tuations whose central part shows crystal-like atomic arrange- 
ment. Those heterophase fluctuations that exceed a critical 
size (determined by the interplay of the interfacial and volu- 
metric contributions to the cluster free energy) have a good 
chance to reach macroscopic dimensions, while the smaller 
ones decay with a high probability. The description of such 
fluctuations is problematic even in single component sys- 
tems. One of the main difficulties is that the typical size of 
the critical fluctuations that form on the human time scale 
CI B Et is comparable to the thickness of the crystal- 
liquid interface, which in turn extends to several molecular 
layers (6J, Q, lU ■ Therefore, the droplet model of the classical 
nucleation theory (CNT) that relies on a sharp interface and 
bulk crystal properties is expected to be inappropriate for de- 
scribing such fluctuations. Field theoretic models, that predict 
a diffuse interface, offer a natural way to handle such prob- 
lems (3- For example, in a recent paper [ilOj, the phase field 
theory (PFT) has been shown to describe such fluctuations 
quantitatively. 

Due to extensive computer simulations done recently, the 
hard sphere (HS) liquid is probably the best known system that 
shows crystal nucleation. Its thermodynamic properties can 
be obtained accurately by integrating the equations of state 
of the crystalline and liquid phases evaluated from molecu- 



lar dynamics (MD) | llj . The interface density profiles are 
known for the (111) and (100) interfaces 1 8]. The free energy 
of the (111), (110), and (100) interfaces has been evaluated 
with a high accuracy 1 12]. Furthermore, the height of the nu- 
cleation barrier has been determined by Monte Carlo (MC) 
simulations 1 13]. This offers a unique possibility to test the 
abilities of various cluster models. For example, it has been 
established that the droplet model of the classical nucleation 
theory (that relies on the equilibrium value of the interface 
free energy) seriously underestimates the height of the nucle- 
ation barrier, and the agreement could only be restored via 
assuming that the interface free energy increases with super- 
saturation |5|, 1 1 311 . Such behavior has been recovered by a 
single-order-parameter density functional approach based on 
a Ginzburg-Landau free energy consistent with face centered 
cubic (fee) crystal symmetries 11411 . Remarkably, however, 
the interface profiles from models deduced for fee symme- 
tries | l^ frxflrill do not fit well to the simulation results. It is, 
therefore, desirable to refine the field theoretic models so that 
details of the interface profiles are reproduced together with 
the height of the nucleation barrier. A possible starting point 
for such study is the phase field theory, that has been used suc- 
cessfully to describe many aspects of crystallization 
including nucleation 1 10]. 

Besides these, further interest in describing crystal nucle- 
ation in the hard sphere system is generated by recent exper- 
iments on colloidal suspensions that mimic closely the hard- 
sphere behavior 1 19]. 

Herein, the hard sphere system is used to test the perfor- 
mance of the phase field theory. The model parameters of the 
phase field theory are fixed at the solid-liquid equilibrium so 
that the known interface thickness and free energy data are re- 
covered. Then, the nucleation barrier height is predicted with- 
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out adjustable parameters, which is compared to exact results 
from Monte Carlo simulations. It will be shown that the phase 
field theory leads to a considerable improvement relative to 
the classical droplet model. Having determined the interface 
density profiles, we evaluate the Tolman length, a quantity re- 
lated to the curvature dependence of the interface free energy. 
We find that in agreement with computer simulations, the Tol- 
man length is positive for small clusters and decreases towards 
a negative value with increasing cluster size. 



II. NUCLEATION THEORY 
A. Phase field theory of nucleation 

In the present work, we apply two fields to describe the lo- 
cal state of matter: A non-conserved structural field m called 
phase field, and the conserved volume fraction field <fi, which 
is related to the density field as <f> = f c 3 p, where a is the HS 
diameter, and p is the number density. For historical reasons, 
to is defined so that it is 1 for the liquid and for the solid. 
Then, 1 — to can be regarded as a local degree of crystallinity, 
and can be viewed as the Fourier amplitude of the dominant 
density wave in the Fourier expansion of the time averaged 
number density in the crystal. 

In line with the formulation for binary systems fTolfTslEol 
l2lll . the Helmholtz free energy of the inhomogeneous hard 
sphere system is assumed to have the form 



F 



dr 



bT 



(Vto) 2 + /(to,0) , 



(1) 



where the local Helmholtz free energy density is given as 

/(m,0 = wTg{m) + [1 - p(m)]fs(4>) + p(m)f L (<f>), (2) 

the g "double well" and p "interpolation" functions have the 
usual form g(m) — jm 2 (l — to) 2 and p(m) = m 3 (10 — 
15to + 6m 2 ), emerging from the thermodynamically con- 
sistent formulation of the phase field theory 1 20, 22], while 
fs(4>) and jx(</>) are the free energy densities of the bulk 
crystal and liquid phases. The square gradient (SG) term in 
the integrand is responsible for the diffuse interface. Being in 
unstable equilibrium, the critical fluctuation (the nucleus) can 
be found as an extremum of this free energy functional. Then 
the fields have to obey the appropriate Euler-Lagrange (EL) 
equations, which in the case of such local functional take the 
form 



5x dx dVx 



(3) 



where ^£ is the first functional derivative of the free energy 

with respect to the field x(= m or <P)> an d V> i s the total local 
free energy density [the integrand of Eq. (1)]. These equations 
have to be solved under the boundary conditions of having un- 
perturbed liquid far from the fluctuation, and due to symmetry 
reasons, zero field gradients at the center of the fluctuations. 
Note, that since the volume fraction is a conserved field, we 



have to add A J drip — — p.^ J drp to the right hand side of 
Eq. (1), when searching for the extremum of the free energy, 
where the Lagrange multiplicator A turns out to be propor- 
tional to the negative of the chemical potential p^ of the un- 
perturbed liquid. This is then the Legendre transformation that 
yields the grand potential, Q = F — PooN, i.e., we seek the 
extremum of the grand potential functional. Assuming spher- 
ical symmetry, the EL equations boil down to the following 



6 Ji = pL-bT{m» + 2 -m'} = (4) 
dm dm r 



and 



sn dip sv> 



(r — > 00) = 0, 



(5) 



where r — > 00 means that the term has to be evaluated far from 
the fluctuation. [In this work ' stands for differentiating with 
respect to the argument where argument is shown; otherwise 
(as above), it denotes differentiation with respect to the radial 
distance r from the center of the fluctuation.] Note that Eq. (5) 
is an implicit equation for the volume fraction = (f>(m), 







and thus 
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= wTg'(m)+p'(m)[f L ((t>) - / s (</>)] 



(6) 



(7) 



is a function of to, i.e., Eq. (4) is an ordinary differential equa- 
tion for m. 

The properties of the critical fluctuation can be determined 
by solving Eq. (4) under the boundary conditions ml — ► for 
r — > 0, and m — ► 1 and <fi — > </>oo for r — ► 00 that corre- 
spond to the initial non-equilibrium liquid state. In this work, 
Eq. (4) has been solved numerically, using a variable 4th/5th 
order Runge-Kutta method. Since the boundary conditions fix 
to and to' at different positions, the phase field value at the 
center, m(r — > 0), for which the far-field condition is satis- 
fied has been found iteratively. 

Having determined m(r), the volume fraction/density pro- 
file can be obtained by inserting it into the numerically in- 
verted version of Eq. (6). The free energy of the critical fluc- 
tuation, in turn, can be calculated by inserting the solution 
m(r) into 



W* 



drAirr 2 Au>[m, (p(m)], 



(8) 



where Aw = ui — lo^, and the grand potential density is u> — 
ip — PocP, while LUoc — ipoo — /iooPoo- Here, subscript 00 
denotes properties of the initial liquid. 

The interface free energy of small fluctuations is expected 
to depend on size or curvature due the reduction of the average 
number of solid neighbors of surface molecules relative to the 
planar interface. The analogous phenomenon in small liquid 
droplets has been studied extensively 12311 . In the widely ac- 
knowledged thermodynamic theory of Tolman ll24Tl . the size 
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dependence of the surface tension is given in terms of the Tol- 
man length 1 25 ] St — R e — Rp ( me distance of the equimolar 
surface Il26ll from the surface of tension |27]) as 



' 1 + 2S T /R P ' 

where is the surface tension for planar geometry. Al- 
though a rigorous derivation of these notions is unavailable 
for crystallites, a quantity analogous to 5j~ has recently been 
evaluated from computer simulations 1 2 ] . It decreases with in- 
creasing size of the fluctuations. It is of considerable interest 
to see whether the phase field theory is able to reproduce this 
feature. 

We determine the Tolman length here via three routes: 

(i) Using the radius of the equimolar surface, R e = 
{f™ drAnr 2 ^ - ( /» oo ]/[(4 7 r/3)(0 o - 0oo)]} 1/3 , where O is 
the value at the center of the fluctuation, and the expres- 
sion R p = [3W*/{2nAuj max )} 1 / 3 O for the radius of 
the surface of tension. Here Auj max = max{fL(<t>oo) + 
(<9/l / 0(f)) ^ (0— 0oo ) -fs (4>) } is the maximum driving force 
of crystallization, realized by the (compressed) solid that has 
the same chemical potential as the initial liquid Il 29ll . 

(ii) On the basis of Eq. (9), Sr, e ff — \Rp{loo/l — 1), 
where the ratio of the curved and planar interfaces is related 
to the ration of the non-classical and classical cluster free en- 
ergies as 7 / 7oo = {W*/W^ NT f^ 0. 

(iii) On the basis of Eq. (9), but approximating the radius 
of the surface of tension with that of the equimolar surface, 
Sx.eff = \Re{loo/l ~ 1)' as done in the evaluation from 
MD results Q]. 

Note that Eq. (9) has been derived from the hypothesis of 
a size-independent Tolman length, a condition that is not sat- 
isfied in the known cases 1 30 ] , thus routes (ii) and (iii) are 
expected to be less accurate. 

Provided that /l(</>) and fs(<f>) are known, the phase field 
theory contains only two parameters, the coefficient of the 
square-gradient term b and the free energy scale w. It is worth 
recalling in this respect that these quantities can be related to 
the free energy and thickness of the equilibrium planar inter- 
faces [[221. 

7oo = (26T) 1 / 2 / d£/E,#0] 1/a , (10) 
Jo 

and 

dio-90 = (^) 1/2 / dtm<f>(or 1/2 , (id 

1 J0.1 

where dio-90 is the 10% - 90% thickness of the interface, the 
distance on which m varies between 0.1 and 0.9. Thus, if 
and dio-90 are known, these relationships can be used to fix 
b and w in equilibrium. This then allows the calculation of 
the free energy of the critical fluctuations W* and the Tolman 
length without adjustable parameters. 

B. Classical theory 

For comparison with our non-classical approach, we cal- 
culate the height of the nucleation barrier using the sharp in- 



terface droplet model of the classical nucleation theory [31]. 
In this approach, the free energy of heterophase fluctuations 
of radius R is given as Wcnt = — (4:^/3) R 3 Au> max + 
47ri? 2 7 oc . Then, the free energy of the critical fluctua- 
tion of radius R% NT — 2^^/ Auj max reads as W£ NT = 
(16^/3) 7 ^A w - 2 a;c . 



C. Application to hard spheres 

The stable crystalline phase of the hard sphere system has 
the face centered cubic (fee) structure. The nuclei, in turn, 
contain fee (ABC) and hexagonal close packed (hep) stacking 
(ABA) with equal probability |5], indicating that the nucle- 
ation barrier for them is rather similar. Therefore, we address 
here only fee nucleation. 

The Gibbs free energy difference between the fee solid 
and liquid (the driving force of crystallization) has been ob- 
tained by integrating the virtually accurate equations of state 
by Hall 1 1 1 ] (that fit to the molecular dynamics results), start- 
ing from the volume fractions of the coexisting solid and 
liquid 0f cc ,cocx = 0.546 and 0l,cocx = 0.494 at pressure 
p = 11.81fcT/(7 3 . Polynomials were then fitted to the chem- 
ical potentials of the solid and liquid phases, whose relative 
deviation from the accurate chemical potentials are less then 
10~ 4 in the density range of interest. These polynomials were 
then integrated analytically to obtain the free energy densities 
of the two bulk phases. 

Since the structural order parameter 1 — m can be inter- 
preted as the Fourier amplitude of the dominant density waves 
(with the first neighbor reciprocal lattice vectors as the wave 
number), we identify the phase field profile as the cross- 
interfacial variation of the height of the density peaks. This is 
an approximation, since other Fourier components also con- 
tribute. However, as pointed out by Shen and Oxtoby fl5ll . 
so far as the density of solid is reasonably approximated by 
Gaussians centered around the fee lattice sites, the amplitudes 
of all other components are proportional to the amplitude of 
the dominant wave, i.e., to a first approximation the height of 
the density peaks is proportional to the amplitude of the dom- 
inant Fourier component. 

The interfacial density profiles of the equilibrium hard- 
sphere crystal-liquid interfaces of orientations (100) and (111) 
have been studied by Davidchack and Laird using molecular 
dynamics simulations |8]. We have received new high resolu- 
tion data from them for the (111), (100), and (110) interfaces. 
In the latter cases, the bin size for the density profile was 1/32 
times the layer width, while the parameters for Gaussian fil- 
tering were N — 64 and e = 28.1 |32]. (For definitions see 
10].) The new density distributions are fully consistent with 
those in [8]. From the interfacial variation of the peak ampli- 
tudes, we obtained 5.94cr, 4.75ct, and 5.54ct for the 10%-90% 
interface thickness for the (111), (110), and (100) directions, 
respectively. 

The same authors evaluated the free energy of the fee crys- 
tal - liquid interface for these orientations 11211 : 7oo<r 2 jkT = 
0.58 ± 0.01, 0.62 ± 0.01, and 0.64 ± 0.01 for the (111), 
(100) and (110) interfaces. Considering that due to the 



FIG. 1: Free energy surface for the hard sphere system evaluated 
with free energy scale w corresponding to the (111) interface. Note, 
that the minima at (0 = 0.494, m = 1) and (<f> = 0.546, m = 
0) correspond to the equilibrium liquid and solid (fee) phases. The 
calculation has been performed for a colloidal suspension of hard 
sphere diameter a = 890 nm, a particle size comparable to that in 
many experiments 1 19]. The small value of the free energy density 
follows from this large molecular diameter. 

small anisotropy the equilibrium shape is expected to be 
nearly spherical, we also performed calculations with an av- 
erage of the interface free energy over the known orientations 
7oo o- 2 /fcr = 0.6133. 

Having this information available, the model parameters 
b and w were determined using Eqs. (10) and (11) via 
Newton-Raphson iteration. The Helmholtz free energy den- 
sity [Eq. (2)] obtained with w evaluated using the properties 
of the (111) interface is shown in Fig. 1. Qualitatively similar 
free energy surfaces were obtained for the other orientations. 

Using these data, the PFT and CNT predictions for the nu- 
cleation barrier can be made free of adjustable parameters. 



III. RESULTS AND DISCUSSION 

A. Density profiles 

Fixing the model parameters b and w as described above, 
the phase field profile for the equilibrium solid-liquid interface 
can be determined via the equation 1 29] 

bT P m 

<m) = (—) 1 / 2 ^/KX0r 1/2 , d2) 

z J0.5 

where z is the distance perpendicular to the plane m = 0.5. 
The coarse grained density profile that corresponds to the fil- 
tered density by Davidchack and Laird [8] can be obtained 
by inserting the solution m(z) into Eq. (6) after inverting it 
numerically. 

The predicted structural order parameter (1 — m) and nor- 
malized coarse grained density (p - p m i n )/(Pmax ~ Pmin) = 



FIG. 2: Reduced interfacial profiles for the (111) interface. (Solid 
lines: predictions of the phase field theory; heavy solid line - struc- 
tural order parameter, 1 — m; light solid line - density. Dashed lines: 
Results from molecular dynamics |8]; heavy dashed line - height 
of singlet density peaks; light dashed line - filtered density profile. 
Note that the heavy lines should be compared to each other, as the 
light ones.) 
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FIG. 3: Reduced interfacial profiles for the (100) interface. Notation 
as for Fig. 2. 

(4> - 4> mm )/ (<pmax ~ 0mm) profiles are compared with then- 
counterparts from molecular dynamics simulations in Figs. 2 
to 4. For the sake of comparison, the simulation results are 
also presented in a normalized form, (X — X min ) / (X max — 
Xmin), where X is either the density peak height or the fil- 
tered density, while X mm and X max are the minimum and 
maximum values of these quantities. 

For the (111) and (100) interfaces, and to a smaller extent 
for the (110) interface, the structural order parameter pro- 
files are in a good agreement with the cross-interfacial vari- 
ation of the density peak height from MD. Note that while 
d 10 -9o has been fitted, the agreement depends also on the 
shape/symmetry of the predicted profile. It appears, that the 
nearly symmetric structural order parameter profile of the PFT 
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FIG. 4: Reduced interfacial profiles for the (110) interface. Notation 
as for Fig. 2. 



approximates reasonably well the behavior seen in the simu- 
lations. Remarkably, the density functional theories by Shen 
and Oxtoby fTslflcjll . and Granasy and Pusztai 1 14] derived for 
the fee structure predict significantly more asymmetric struc- 
tural order parameter profiles for the crystal-liquid interfaces, 
which are sharper on the crystal side and have a long tail on 
the liquid side. Therefore, their match to the simulated pro- 
files is apparently less satisfactory. It is an intriguing question 
why the quartic form of g(m), that is consistent with the sym- 
metries of the base centered cubic (bec) structure y^] , leads 
to a better fit to the simulation profile then those obtained by 
using free energy functionals consistent with the fee symme- 
tries. A distinct possibility is the presence of a bcc-like layer at 
the interface as reported in the case of the Lennard- Jones sys- 
tern HQ. This question, however, could only be addressed 
in a full density functional theory (such as that by Shen and 
Oxtoby 1 34]) that incorporates the possibility for the Bain dis- 
tortion that connects the bec and fee structures. 

The coarse grained density profiles can now be predicted in 
the PFT without adjustable parameters. They appear to be in a 
good agreement with the simulation results for the (111) and 
(100) interfaces (Figs. 2 and 3), however, a shift of about half 
a a can be seen between them for the (110) direction, the PFT 
result lying closer to the crystal (Fig. 4). 

Summarizing, in the PFT, we have a free energy surface and 
a SG coefficient that reproduce exactly the thermodynamics 
of the bulk phases, the interface free energy, and quite reason- 
ably both the structural order parameter profile and the den- 
sity profile. Thus, we may attempt to calculate the height of 
the nucleation barrier with some confidence. 



B. Nucleation 

The reduced nucleation barrier heights calculated with b 
and w that fit to the properties of the (111) and (110) interfaces 
and to the average interface properties are shown in Fig. 5 for 
the phase field theory. For comparison, the predictions of the 
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FIG. 5: Reduced nucleation barrier height vs. volume fraction of the 
initial liquid. (Solid lines: phase field theory - PFT; dashed/dotted 
lines: classical droplet model - CNT.) The upper an lower curves 
were calculated using the physical properties of the (111) and (110) 
interfaces, respectively. The prediction for the (100) interface (not 
shown here) falls slightly below the results obtained with the aver- 
age interface properties (heavy lines). For comparison, the results of 
Monte Carlo simulations [ 13] are also shown (triangles). 



classical nucleation theory and the exact results from Monte 
Carlo simulations 11311 are also presented. The W* vs. initial 
liquid volume fraction curves predicted by the PFT using the 
properties of the (111) and (110) interfaces envelope the MC 
simulations, while the predictions with the average interface 
properties fall close to the MC results. The slope of the W* 
vs. ifioc curve appears to be somewhat larger for the PFT pre- 
dictions than for the MC simulations. A possible explanation 
could be the density dependence of the parameters b and w 
(3- However, the more complex Euler-Lagrange equations 
associated with such problem will be investigated elsewhere. 

It is of interest to determine the magnitude of the error asso- 
ciated with these results. While the relevant thermodynamic 
data are virtually exact, uncertainties of (±0.01fcT/cr 2 ) and 
(±5%) are associated with the interface free energy and the 
interface thickness. We find that for <f> = 0.5207 the barrier 
height is W* jkT — 40.63, while its uncertainties originat- 
ing from those of the interface free energy and the interface 
thickness amount to ±1.91 and ±0.23, respectively, yield- 
ing ±1.92 when assuming Gaussian error propagation. The 
magnitude of these errors varies approximately proportionally 
with the barrier height. These errors do not influence the va- 
lidity of our conclusions. 

The radial order parameter and density profiles correspond- 
ing to the volume fractions used in MC simulations (0.5207, 
0.5277, and 0.5343) are shown in Fig. 6. The interfacial pro- 
files for = 0.5207 are in a reasonable agreement with 
the size of the respective critical fluctuation from MC simu- 
lation (cf. our Fig. 6 and Fig. 3 of Ref. |5|]). Note that the 
bulk crystalline properties have not been established even at 
the centers of the critical fluctuations. Accordingly, the clas- 
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FIG. 6: Radial structural order parameter (solid lines) profiles and re- 
duced density (dashed lines) profiles at initial liquid volume fractions 
(j>oo = 0.5207, 0.5277, and 0.5343, calculated with b and w corre- 
sponding to the average interface properties (7ooO" 2 /fcT = 0.6133 
and dio-go/fT = 5.410). Here S (= 0.5796, 0.5875, and 0.5948, 
respectively) is the volume fraction of the solid that is in mechani- 
cal equilibrium with the initial liquid (the crystal that has the same 
pressure as the liquid of volume fraction 4>oa). 



sical droplet model significantly underestimates the height of 
the nucleation barrier for all orientations. Since the nucle- 
ation rate is proportional to exp(—W*/kT), these differences 
amount to orders of magnitude. For example, the calculation 
made with the average interface free energy overestimates the 
nucleation rate by three to five orders of magnitude. At the 
volume fraction (p^ = 0.5781, the PFT and CNT predictions 
intersect each other. In both approaches, the height of the 
nucleation barrier decreases monotonically with the supersat- 
uration. 

The three known orientations can in principle be used to 
approximate the anisotropy of the interface free energy 1351 
13611 . that in turn defines the equilibrium shape of minimum 
interface free energy for given volume. Since the nuclei are 
expected to minimize their free energy, their average shape 
can probably be reasonably approximated by the equilibrium 
form. Work is underway to determine the nucleation barrier 
in such case. 

Finally, we draw attention to the fact, that the MC data for 
used here as reference, were obtained directly 
via "umbrella sampling", therefore, possible uncertainties as- 
sociated with the nucleation prefactor yl do not plague them. 



C. Tolman length 

The Tolman lengths calculated for the critical fluctuations 
along the three routes specified above are shown as a func- 
tion of the initial volume fraction in Fig. 7. Similarly to 
the field theoretic results for vapor- liquid nucleation |3jj|, for 
small sizes the Tolman length is positive and comparable to 
the molecular diameter, while 6t decreases with increasing 



1.5 




0.5 0.55 0.6 0.65 



4> 

oo 

FIG. 7: Tolman length as a function of initial volume fraction calcu- 
lated along routes (i)-(iii). 



size towards a negative limiting value for the planar inter- 
face. Although this behavior is general for the Tolman lengths 
from all three routes, routes (ii) and (iii) yield considerably 
smaller values than the more accurate route (i). Accordingly, 
the change of sign happens at different volume fractions: At 
4>oo = 0.5266 for route (i), while at = 0.5781 for routes 
(ii) and (iii). The differences originate from the fact that the 
assumption St = const., routes (ii) and (iii) rely on, is not 
satisfied for small clusters. In turn, in the large particle limit 
this assumption becomes valid. Indeed, the Tolman lengths 
from all three routes converge to the same negative value 
($T,eq/G = —0.33) in the planar limit. 



IV. SUMMARY 

We presented a phase field theory of homogeneous crystal 
nucleation in hard sphere liquid, that describes local state of 
matter via a structural order parameter field and the density 
field. 

A rigorous test of the model is performed using the virtually 
exact thermodynamic properties of the bulk phases. The two 
model parameters are chosen so that the free energy and the 
10% — 90% thickness for the cross-interfacial variation of the 
density peak height known from molecular dynamics simula- 
tions are recovered. For the (111) and (100) interfaces, and to 
a lesser extent for the (110) interface, the resulting structural 
order parameter and density profiles are in a significantly bet- 
ter agreement with the computer simulations, than seen for 
previous continuum models. Without adjustable parameters, 
the phase field theory predicts the height of the nucleation bar- 
rier with a reasonable accuracy. This represents a consider- 
able improvement over the sharp interface droplet model of 
the classical nucleation theory, which is known to overesti- 
mate the nucleation rate by several orders of magnitude. 

The cross-interfacial density profiles are used to calculate 
the Tolman length for critical fluctuations, which is positive 
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for small sizes, and tends to a negative value in the large clus- 
ter limit. 
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